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Abstract 

We calculate the Lyapunov exponents in a classical molecular dynamics 
framework. The system is composed of few hundreds particles interacting 
either through Yukawa (Nuclear) or Slater-Kirkwood (Atomic) forces. The 
forces are chosen to give an Equation of State that resembles the nuclear 
and the atomic 4 He Equation Of State respectively near the critical point 
for liquid-gas phase transition. We find the largest fluctuations for an initial 
"critical temperature". The largest Lyapunov exponents A are always positive 
and can be very well fitted near this "critical temperature" with a functional 
form A oc |T — T c |~ w , where the exponent u = 0.15 is independent of the 
system and mass number. At smaller temperatures we find that A oc T °- 4498 ) 
a universal behavior characteristic of an order to chaos transition. 

PACS numbers: 5.45b, 5.70Fh, 24.60Lz 
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The phenomena of phase transitions have been always a subject of great interest for many 
(generations of) physicists. Very recently there has been a large interest in studying what 
happens when the system is not composed of 10 23 particles and confined in a box but, on 
the opposite, is composed of few hundred particles that are not confined. This interest has 
born first in heavy ion collisions where one would like to study the Nuclear Equation of 
State (NEOS) [|lj. In the nuclear case it is not possible to study the infinite number of 
particles case (apart in stars with the obvious difficulties), so one feasible way is to perform 
proton-nucleus or nucleus-nucleus collisions. In this manner it is possible to heat up and to 
explore different densities of the nuclei. Of course the problem is complicated apart from the 
fact of having a small number of particles, also from the presence of Coulomb force, angular 
momentum and from the approximate knowledge of the "thermodynamical parameters", i.e. 
temperature, density and pressure @-[|]. Similar problematic can be found in the study of 
metallic clusters and/or fullerenes ||. One of the first questions that arises is: does it make 
sense to speak about phase transitions for a system made of 100-400 particles? To answer 
this question the authors of refs. and |5| have solved the exact classical many body 
problem (Classical Molecular Dynamics-CMD) numerically for particles interacting through 
two body forces. In particular preparing 100 particles in the ground state and giving to the 
particles an excitation energy (or temperature T) the following results are found: 

1) for a given initial temperature T c the system undergoes multifragmentation [|]||[7|]. 
The resulting mass distributions display a power law Y(A) oc A~ T with r = 2.23, which 
is exactly the value predicted in the Fisher droplet model for a system near the liquid-gas 
phase transition point ||. Other critical exponents have also been estimated within the 
same model 0:7 = 0.86, (3 = 0.33 [|]. Similar estimates but for experimental data on 
Au+C collisions at 1 GeV/nucleon Q give also reasonable values of all the exponents; 

2) at the same temperature, fluctuations in the mass distribution are maximum. This is 
revealed through the study on an event by event analysis of intermittency signal, factorial 
moments and Campi's plots [p| JTTJ| , [TT| ] ; 

3) the "critical temperature" follows the scaling law T c /\Eb\ = constant ||. Where Eb 
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is the binding energy of the system (10.5 MeV for nuclei, 9.5 eV for C60 fullerenes and 50.5 
°K for A He molecules); 

4) at very small temperatures or, equivalently, small excitation energies, the events are 
typical evaporation like events, i.e. with a final mass distribution composed of a big fragment 
and many small ones : monomers, dimers etc. 

For some initial conditions the system displays large fluctuations, thus we expect that 
other indicators of fluctuations -the Lyapunov exponents - should be positive. The values 
of the Lyapunov exponents for systems undergoing a phase transition are not known (yet). 
In particular the relationship between thermodynamical and dynamical properties have not 
been exhaustively explored. Ours is the first estimate of the Lyapunov exponents for sys- 
tems exhibiting a critical behavior (possibly because of a liquid-gas phase transition). In 
a previous exploratory work, but for a two dimensional system and in the mean-field ap- 



proximation, the largest Lyapunov exponents were calculated in the spinodal region |T2 
It is well known, however, that the mean-field approximation gives a rough estimate of the 
critical exponents ||, thus it is interesting to estimate their values in CMD. 

We can summarize the main results of this work as follows: 
independently of the studied system, i.e. nuclei, helium molecules, the largest Lyapunov 
exponents satisfy, similarly to the Landau theory of phase transitions, the relation 

A oc |T-T c p, T~T C (1) 

where u = 0.15 ±0.04. 

At smaller temperature the Lyapunov exponents scale as 

A oc rra, (2) 



with 5 = 4.669... the Feigenbaum constant [OA which indicates a typical transition from 



order to chaos . 

Our studies are based on direct numerical simulation of an expanding system in classical 
molecular dynamics. In particular we have studied a "nuclear" system where "neutrons" 



and "protons" interact through Yukawa potentials. Details of the forces can be found in 
refs. [||,[nj]. The parameters entering the potentials have been fitted in such a way to have 
a ground state density of 0.16 fm~ 3 , a binding energy of -16 MeV/nucleon for an infinite 
system. For finite systems the binding energy is -10.5 MeV/nucleon without Coulomb. In 
this work we will discuss the parameter set that gives a soft EOS, i.e. compressibility K=200 
MeV |16| and the Coulomb interaction will be neglected. We notice that the use of the set 



corresponding to a stiff EOS, K=380 MeV, gives the same results. 

We have also studied He atoms interacting through the potential [17 

v(r) = 5.67 ■ 10 6 e - 2L5(r/<T) - 1.08 (^) 6 , (3) 

where a = 4.64 Angstrom and v(r) is in Kelvins. 

In order to calculate the largest Lyapunov exponents [IE] we first define the following 
metric in phase-space 



d(t) 



N 



\ £[« (r'i(t) - r' 2 (t)y + (3 (p^t)) - p' 2 (t)n , (4) 
\ i=i 



where the sum runs over all the N particles of the system, the subscript 1 and 2 refers to 
two events that at time t=0 differ of an infinitesimal quantity d(0) = 10 -6 or less. The r' 
and p' are scaled positions and momenta. In our case 

r'(t) = r(t)/R rms , (5) 

p'(t) = p(t)/P , (6) 

where R rm s is the root mean square radius and P = \ / 2mT is an average momentum, being 
m the nucleon (or the atom) mass and T the initial temperature (see below) . Normally the 
Lyapunov exponents are calculated for systems bound in phase-space. This is not always 
the case in our simulations since for high excitation energy the phase-space explodes. In 
order to be certain of the soundness of our results we calculated the metric for three cases: 
1) a=l and /3=1; 2) a—1 and (3=0; 3) a=0 and (3=1. The results are independent on the 
choice of the metric, cases 1-3, as it should be. 



The Lyapunov exponents A are obtained from the relation 

d(t) = d(0) e xt . (7) 

In our numerical simulations we prepared a system of 100 particles in their ground state. 
Then we distributed the momenta of the particles according to a Maxwell-Boltzmann distri- 
bution at temperature T and let the system evolve in time by solving the classical equations 
of motion. At each temperature 100 events were generated. For each event (test event) ten 
other events were generated, each event differing from the test event of d(0). The exponents 
were obtained by averaging over all the events. In order to get A numerically, we calculated 
the ratio d(t)/d(0) and fitted its exponential growth. As a check we estimated also 

d(t)/d{t) = A . (8) 

In figure 1 we plot typical evolutions of d(t)/d(0) at three temperatures T=2, 5 and 20 MeV 
for the nuclear system. We see that the distance increases exponentially with time and can 
be very well fitted with a straight line (in a semilog plot) whose slope is just the Lyapunov 
exponent. The same results are found by using eq.(8). In general we get Lyapunov exponents 
which are always positive at all temperatures and have a maximum value at T ~ 4.5 MeV. 
A similar behavior is found for the atomic case with a maximum at T ~ 21.8 °K. In order 
to compare such different systems like nuclei or atomic He, we scale temperatures and times 
with typical values of the different systems. Usually one scales with the values obtained at 
the critical point. Since we do not know these values a priori, we scale the temperature with 
the absolute value of the binding energy and the (inverse) times with a typical value Ao 

A = p\E b \/m/R . (9) 

Recall that in order to derive the EOS || one needs to know the hard core radius and 
the depth of the two body potential. These quantities are proportional to the average 
equilibrium distance between particles Rq = R/A^, - where R is the radius of the system 
and A the mass number - and to the binding energy. In our case Ao = 0.155 c/ fm for 
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the nuclear case and Ao = 8.82 • 10~ 5 ac/ciQ for the atomic one, being a the fine structure 
constant, c the speed of light and ao the Bohr radius. In figure 2 we plot the Lyapunov 
exponents at each initial temperature (in units of their typical values) for the nuclear (circle 
symbols) and atomic cases (squares). Note that the maximum value of A is obtained at the 
initial scaled temperature 

T c /\E b \ ~ 0.43 (10) 

for both systems. For such temperatures the systems undergo multifragmentation and a 
power law in the mass distribution and factorial moments is found [§,|J. Thus our result 
confirms that in this temperature region fluctuations are largest. Inspired by the Landau 
theory of phase transition, we parametrized the exponents according to eq.(l). The fit is 
also displayed in figure 2 (full and dashed curves). The same value u = 0.15 ±0.04 was 
used in both cases. Note the good agreement with the scaled A at all temperatures but 
the lowest ones. For such small temperatures there is not any multifragmentation of the 
system, indeed the mass distribution reveals a typical cases of evaporation. Using the same 
language as in the theory of phase transitions we could say that the small temperature cases 
explore densities and pressures outside the spinodal region. Note that the absolute values 
of the scaled exponents differ of less then 30% for the two cases. Correspondingly the two 
fits differs only for a multiplicative constant C. The reason for such a small discrepancy 
is due to the fact that in the nuclear case we have two fluids, neutrons and protons. We 
will discuss this point in more detail in a following publication |Tj|. We also stress that 
these results are independent on the chosen metric, i.e. in the values of a and j3 in eq.(4). 
We have also tested the results by changing the mass of the fragmenting system. The 
Lyapunov exponents remained the same for masses ranging from 50 to 400 particles |jTH . 
The uncertainties reported in figure 2 are of the order of ~ 10%. 
A scaling law of Lyapunov exponents of the kind 

A oc (A - A*) ln2 ' lnS = (A - A*) - 4498 , (11) 
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where A is a control parameter and S = 4.669... the Feigenbaum constant [13|, is typical 



of order to chaos transitions. The value A* is the critical value which indicates the onset 



of chaos. The expression (11) was initially found in the logistic map |TJj], but various 
experiments have confirmed its general validity ||15|| . In our case the largest Lyapunov 
exponents A are positive for all finite temperatures, i.e. the dynamics is always chaotic, but 
they tend to zero as the temperature T goes to zero. Actually for T=0 the systems are 
frozen in their ground state and A vanishes. Thus chaos starts at T > 0. 

In figure 3 we plot a magnification of figure 2 for very small temperatures. The full curve 
is the expression (11) multiplied by a constant fitted on the numerical points. In this case 
the control parameter is the scaled temperature and A* corresponds to T=0. The agreement 
is really impressing for both systems up to T/\Ef,\ ~ 0.1. 

We can try to understand this behavior by recasting our microscopic dynamics in terms 
of a phenomenological macroscopic model. At very low temperatures one gets essentially 
evaporation events. Considering a discrete map, at time n the system evaporates z n particles 
and it will continue evaporating qz n particles until the excitation energy is zero. The number 
of evaporated particles at time n+1 is thus z n+ \ — (1 + q)z n . If z is the final number of 
particles and we assume that the number of evaporated particles decreases linearly with time 
because the excitation energy is decreasing, we obtain z n+ \ = (1 + q)z n {l — z n /z ). This is 
nothing else that the logistic map if 1 + q — A and x n = z n /z . 

Thus we have two different mechanisms at play. The first one for small temperatures 
gives a transition from order (the ground state) to chaos and has a dynamical origin. The 
second mechanism, for reduced temperatures larger then 0.1 has a thermodynamical origin 
appropriate for a second order phase transition. Loosely speaking we have given evidence 
for "critical chaos" in the latter case. 

In conclusion, in this work we have calculated the largest Lyapunov exponents as a 
function of the initial temperature for an expanding system composed by 100 particles in 
the framework of classical molecular dynamics. We have shown that these exponents are 
always positive and have their largest value at a temperature of ~ 4.5 MeV for the nuclear 



case and ~ 21.8 °K for the atomic one. We have also demonstrated that the A obey 
to universal scaling laws. They fulfill the relation |T — T c |~ w , uo = 0.15 ± 0.04, similarly 
to the Landau theory of phase transitions near the critical point. At the same time for 
smaller temperatures (evaporation events) Lyapunov exponents exhibit a general transition 
from order (T = 0) to chaos (T > 0). We feel that further investigations following the 
ideas presented in this paper may help our understanding of order and disorder in classical 
systems and, after all, in (part of) Nature itself. 
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FIGURES 

FIG. 1. The ratio d(t)/d(0) is plotted as a function of time at three initial temperatures 
T = 2,5, 20 MeV for the nuclear system. The dashed lines are fits whose slope give the typical A 
for these temperatures after averaging over hundreds of events. 

FIG. 2. The scaled largest Lyapunov exponents A/Ao are plotted vs. the scaled initial tem- 
perature T/\Ef,\ for the nuclear (circles) and the atomic (squares) case. The full and dashed lines 
are fits obtained with the functional form C\T — T c \~ w where uj = 0.15 ± 0.04. The parameters of 
the fits are C=0.25 and C=0.3 for the nuclear and atomic case respectively. See text for further 
details. 

FIG. 3. A magnification of figure 2 at very small temperatures. The full curve is a fit with the 
functional form K T 4498 where K = 0.55 is the fitted parameter. See text for further details. 
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